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Abstract 


High  maneuverability  is  one  of  the  major  goals  in  rotorcraft  design.  In  practice, 
this  goal  is  limited  by  unsteady  (dynamic)  stall  near  blade  leading  edges.  Studies 
of  three-dimensional  boundary  layer  separation  on  a  rotating  blade  are  made.  For 
hovering  flight,  the  blade  twist  and  downwash  are  included  in  the  effective  angle  of 
attack.  For  forward  flight,  high  angles  of  attack  are  used  to  simulate  the  most  se¬ 
vere  situation  at  the  retreating  blade.  Because  of  the  disparate  scales  of  the  leading 
edge  radius  and  the  blade  span,  separation  is  found  to  be  quasi  two-dimensional, 
and  local  singular  behaviors  at  separation  are  very  similar  to  the  two-dimensional 
case.  Most  of  the  results  are  obtained  using  an  Eulerian  approach,  but  a  Lagrangian 
formulation  is  used  to  study  the  behavior  near  the  separation  singularity.  Control 
mechanisms  based  on  suction  and  blade  oscillations  are  examined.  It  is  found  that 
oscillations,  with  a  tuned  frequency  and  amplitude,  can  delay  separation.  Leading 
edge  suction/injection  is  also  effective  in  delaying  separation  for  particular  (opti¬ 
mized)  slot  locations. 
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Chapter  1 
Introduction 


1.1  Background 

Dynamic  stall  refers  to  a  series  of  events  on  pitching  airfoils  during  unsteady 
flight.  Stall  is  usually  initiated  by  boundary  layer  separation  near  the  leading  edge, 
and  culminates  in  the  formation  of  a  dynamic  stall  vortex  accompanied  by  a  strong 
interaction  between  the  viscous  inner  flow  and  the  inviscid  outer  flow.  This  can 
result  in  a  severe  loss  of  lift  and  an  abrupt  increase  in  pitching  moment.  During 
changing  flight  conditions,  dynamic  stall  plays  a  decisive  role  in  restricting  airfoil 
performance.  In  experiments,  it  is  usually  observed  that  a  primary  vortex  forms  and 
resides  on  the  upper  surface  of  an  airfoil,  accompanied  by  a  considerable  increase 
in  the  lift.  This  extra  gain  in  lift,  however,  is  not  sustained.  A  secondary  or 
even  tertiary  vortex  is  quickly  seen  beneath  the  primary  one  and  the  interaction 
between  these  vortices  ultimately  causes  the  primary  vortex  to  leave  the  surface 
and  travel  downstream.  As  a  result,  the  flow  held  is  strongly  disturbed  and  the 
airfoil  performance  is  severely  limited. 

Dynamic  stall  is  the  major  impediment  to  enhancing  the  performance  of  rotating 
blades.  Currently,  two  approaches  are  used  to  control  dynamic  stall.  One  approach 
is  involved  with  keeping  the  stall  vortex  on  the  airfoil  surface  so  that  the  benehcial 
aspect  of  dynamic  stall  can  be  retained.  Marginal  separation  theory  may  shed  some 
light  in  this  direction,  but  adequate  theories  or  effective  experimental  methods  are 
not  yet  available.  The  second  approach  is  concerned  with  techniques  for  suppressing 
or  inhibiting  separation  at  high  angles  of  attack.  These  techniques  are  based  on 
delaying  the  vorticity  eruption  that  gives  rise  to  stall. 

Unsteady  separation  control  cannot  be  achieved  in  a  conventional  way.  A  pri¬ 
mary  difficulty  is  that  the  boundary  layer  near  the  leading  edge  evolves  over  a  local 
time  scale  that  is  much  smaller  than  a  global  time  scale  associated  with  the  chord 
length  (airfoil  pitching,  wake  shedding  etc.,  are  all  related  to  the  chord  length).  A 
global  time  scale  [T]  can  be  dehned  as  the  interval  in  which  a  blade  advances  a  half 
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chord.  As  an  specific  example,  if  the  chord  length  c  =  0.5  ft  and  the  free  stream  ve¬ 
locity  Uo  =  350ft/s,  then  the  global  time  size  is  [T]  =  A  local  time  can 

be  estimated  by  [t]  =  ff-,  where  Vg  =  Roe^{c/2),  is  the  radius  of  the  leading  edge,  e  is 
the  thickness  length  ratio  of  an  airfoil  and  Rg  is  a  constant.  For  NACA0012  airfoil, 
e  =  0.12  and  Rg  =  1.1.  Simple  arithmetic  shows  that  the  absolute  value  of  [t]  is  very 
small,  say  around  10“®  second;  it  is  less  than  1/70  of  the  global  scale.  Therefore, 
any  event  inside  the  boundary  layer  is  difficult  to  capture  in  experiments,  especially 
in  a  helicopter  environment  involving  a  variety  of  complex  fluid  phenomena. 

Another  reason  for  the  failure  of  a  conventional  strategy  lies  in  the  fact  that 
boundary  layer  separation  initiates  in  very  small  spatial  scales  above  the  airfoil 
surface.  Even  at  the  instant  prior  to  separation,  measurable  features  on  the  wall 
such  as  surface  pressure  or  wall  shear  do  not  signal  the  onset  of  separation.  Thus 
methods  for  separation  prevention  can  not  be  actuated  by  a  precursor  detection.  In 
order  to  be  effective,  control  should  interfere  at  an  early  time. 

A  closer  look  at  the  physics  of  unsteady  separation  may  be  helpful  before  im¬ 
plementing  any  control  mechanism.  The  development  of  unsteady  separation  near 
the  leading  edge  of  an  airfoil  can  be  characterized  by  two  stages.  In  the  first  stage, 
the  adverse  pressure  gradient  on  the  upper  surface  causes  gradual  deceleration  of 
the  flow  inside  the  boundary  layer  and  strong  vorticity  diffusion  at  the  wall.  In 

dll 

two-dimensional  flow  the  vorticity  u  ~  —  is  initially  negative  everywhere.  The 

dy 

vorticity  flux  across  the  surface  is 


duj 


dx  dt  dx 


VyjUJ 


y=0 


(1.1) 


For  a  solid  surface,  the  vorticity  flux  is  determined  by  the  pressure  gradient,  and 
du) 

<  0  when  the  pressure  gradient  is  positive  (adverse).  The  initially  negative 

vorticity  is  gradually  cancelled  out  by  this  flux  and  a  zero-vorticity  line  finally  ap¬ 
pears.  Development  of  the  boundary  layer  then  comes  to  a  new  stage  in  which  the 
zero-vorticity  line  plays  an  essential  role.  At  this  time,  the  flow  begins  to  recirculate 
in  certain  regions  and  the  shear  stress  is  zero  at  points  on  the  wall.  Separation, 
however,  does  not  occur  at  this  stage.  The  boundary  layer  approximation  is  still 
valid  and  the  displacement  effect  of  the  viscous  layer  on  the  outer  inviscid  flow  is 

At  the  second  stage,  the  zero-vorticity  line  is  elevated  by  convection.  In  the  re¬ 
circulation  region,  convection  dominates  diffusion  and  the  flow  is  essentially  inviscid. 
The  momentum  equation  can  be  approximately  written  as 


du  du 
dt  dx 


0 


(1.2) 


which  is  a  one-dimensional  Burgers  equation.  In  the  inviscid  flow,  the  zero-vorticity 
line  corresponds  to  a  material  line.  Fluid  particles  trapped  on  this  line  have  different 
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speeds.  At  a  finite  time,  particles  on  the  line  collide  at  a  point,  forming  a  singularity. 
A  spike  appears  in  the  displacement  thickness  and  the  boundary  layer  equations 
break  down.  This  is  known  as  the  Van  Dommelen-Shen  process  (see  Van  Dommelen 
&  Shen  1980,1982),  and  is  characterized  by  increased  thinning  along  the  streamwise 
direction  and  thickening  in  the  normal  direction.  The  influence  of  the  viscous  effects 
are  now  much  larger  than  0{Re~^^'^).  A  large  amount  of  vorticity  is  shed  into  the 
outer  flow  over  a  very  narrow  region  in  the  streamwise  direction.  Doligalski,  Smith  & 
Walker  (1994)  conjectured  that  the  vorticity  which  erupts  from  the  boundary  layer 
rolls  up  into  the  primary  stall  vortex,  leading  to  a  series  of  complex  phenomena. 

One  objective  of  the  current  research  is  to  simulate  unsteady  boundary  layer 
separation  in  a  three-dimensional  environment.  As  a  helicopter  blade  usually  has  a 
large  aspect  ratio,  and  near  the  leading  edge  variations  in  the  spanwise  direction  are 
smaller  than  streamwise  variations,  the  three  dimensional  separation  to  a  certain  ex¬ 
tent  is  similar  to  two-dimensional  separation.  The  spanwise  velocity  has  an  influence 
on  separation,  but  the  effect  is  not  large.  However,  the  separation  structure  in  three 
dimensions  is  more  complex.  For  instance,  in  contrast  to  the  two-dimensional  case, 
the  appearance  of  recirculating  flow  in  a  plane  cut  along  any  direction  does  not  nec¬ 
essarily  indicate  separation.  Since  the  vorticity  now  has  two  principal  components 
dvj  dw 

(Jz  ~  —7—  and  (Jx  ~  7:—,  separation  is  initiated  on  the  zero- vorticity  line  associated 
oy  ay 

with  the  intersection  of  the  two  iso-surfaces  Ux  =  0  and  Uz  =  0.  Due  to  blade  twist 
in  flight,  the  angle  of  attack  varies  along  the  spanwise  direction.  The  rotor  trim 
scheme,  the  tip  vortex,  and  the  inboard  trailing  wake  make  the  flow  environment 
highly  unsteady  and  simulation  based  on  a  comprehensive  model  is  challenging.  A 
trailing  wake  model  can  be  constructed  by  extending  the  two-dimensional  model 
of  Zalutsky  (2000)  into  three  dimensions.  In  this  model,  a  wake  is  shed  from  the 
trailing  edge  of  a  thin  airfoil  undergoing  unsteady  motions.  The  principle  parameter 
is  the  strength  of  the  wake  vorticity  7(t),  which  equals  the  tangential  velocity  jump 
at  the  trailing  edge.  7  satishes  an  integral  equation  which  can  be  solved  once  the 
temporal  variation  of  the  angle  of  attack  is  specihed.  Unsteady  loads  can  then  be 
evaluated.  In  three-dimensional  flow,  7  =  'y{z,t)  also  satishes  an  integral  equation 
(Katz  &  Plotkin  2001  ).  For  a  helicopter  with  a  high  aspect-ratio  blade,  the  integral 
equation  can  be  solved  in  a  much  more  efficient  way  by  combining  lifting  line  theory 
and  the  numerical  procedures  developed  by  Zalutsky  (2000).  Related  asymptotic 
methods  can  be  developed  for  rotating  blades  to  hnd  the  leading  order  inhuence  of 
the  unsteady  wake  on  leading  edge  separation. 

For  large  aspect  ratio  wings,  the  separation  process  is  quasi  two-dimensional, 
and  re-examination  of  two-dimensional  results  is  very  helpful.  (More  general  three- 
dimensional  separation  structures  were  considered  by  Van  Dommelen  &  Cowley 
(1990)  Two-dimensional  unsteady  separation  and  control  mechanisms  are  outlined  in 
Chapter  2.  Control  studies  for  oscillating  wings  and  for  leading  edge  suction/injection 
are  examined  in  detail.  Results  of  the  two-dimensional  simulations  are  discussed  in 
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Chapter  3.  It  is  shown  that  there  exists  a  band  of  frequences  and  amplitudes  over 
which  the  separation  time  is  signihcantly  increased.  Suction/injection  calculations 
indicate  that  there  can  be  an  optimal  location  for  the  suction/injection  slots  that 
also  increases  the  separation  time.  Unsteady  three-dimensional  boundary  layer  sep¬ 
aration  is  discussed  in  Chapter  4.  Both  Eulerian  and  Lagrangian  formulations  are 
presented.  Results  of  the  three-dimensional  computations  are  given  in  Chapter  5. 
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Chapter  2 

Two-dimensional  Unsteady 
Separation  and  Control 

2.1  Introduction 

As  noted  earlier,  for  helicopter  blades  with  large  aspect  ratios,  unsteady  sepa¬ 
ration  is  quasi  two-dimensional.  Consequently,  separation  and  control  studies  for 
two-dimensional  flows  provide  useful  guidance  for  prevention  of  three-dimensional 
separation  on  blades. 


2.2  Unsteady  separation  without  control 

For  an  airfoil  at  incidence,  because  of  a  favorable  pressure  gradient,  the  velocity 
hrst  accelerates  along  the  upper  surface.  Leading  edge  curvature,  however,  varies 
in  such  a  way  that  an  adverse  pressure  gradient  is  encountered  after  a  short  dis¬ 
tance.  Within  the  boundary  layer,  fluid  particles  have  a  much  smaller  momentum 
compared  with  the  outer  flow,  but  have  to  overcome  the  same  resistance  due  to  the 
adverse  pressure  gradient.  Beyond  a  critical  incidence,  boundary  layer  separation  is 
expected.  In  unsteady  flows,  flow  reversal  or  zero  wall  shear  is  no  longer  the  dehning 
characteristic  for  separation.  A  zero  vorticity  line,  which  appears  in  the  recirculat¬ 
ing  flow,  signals  the  onset  of  the  Van  Dommelen-Shen  process  (Van  Dommelen  & 
Shen  1980,1982).  In  a  local  zone,  the  boundary  layer  scale  increases  in  the  normal 
direction  and  shrinks  in  the  streamwise  direction.  After  a  hnite  time,  a  singularity  is 
generated  by  the  collision  of  particles  on  the  zero  vorticity  line  and  the  subsequent 
boundary  layer  eruption  ejects  a  large  amount  of  vorticity  into  the  outer  layer.  The 
boundary  layer  approximation  then  breaks  down  and  a  strong  viscous/inviscid  in¬ 
teraction  associated  with  the  next  stage  has  to  be  described  by  another  subset  of 
the  Navier-Stokes  equations.  Walker  (2002). 

Boundary  layer  separation  initiates  at  the  leading  edge.  Analyses  within  this 
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region  are  important  for  control  studies.  In  notes  by  Van  Dyke  (1956),  the  airfoil 
leading  edge  zone  is  approximated  by  a  parabola.  Thin  symmetric  Joukowski  airfoils 
provide  a  simple  example.  They  are  dehned  by 


V 


4e 


(2.1) 


where  X  =  — 1  is  at  the  leading  edge,  X  =  1  is  at  the  trailing  edge,  e  denotes  the 
thickness  ratio  of  the  airfoil,  and  (X,  V)  are  dimensionless  Cartesian  coordinates 
with  respect  to  the  half-chord  length  c/2.  The  equation  of  a  parabola  with  the  same 
vertex  point  and  dimensionless  nose  radius  is 

V  =  v/2r„(X  +  l),  (2.2) 


and  To  =  64e^/27  for  the  Joukowski  airfoil  dehned  in  (2.1).  Local  coordinates  near 
the  nose  may  be  recast,  with  repect  to  r^,  as 


X  =  (X  +  l)/ro,  y  =  Y/to 


(2.3) 


In  the  inviscid  how,  the  surface  velocity  on  the  parabola  at  incidence  is  given  by 
Van  Dyke  (1956). 


M 

A 


X 


XT  1/2 


(2.4) 


a'  is  proportional  to  the  geometric  incidence  a  (see  2.6  below).  For  a  thin  airfoil  at 
incidence  a,  the  velocity  away  from  the  leading  edge  is  (Katz  &  Plotkin  2001) 


L.l  +  ii(l_2A-)±a.y^  (2.5) 

to  hrst  order.  Uo  is  the  uniform  velocity  in  the  free  stream  and  a  is  the  angle  of 
attack.  Matching  with  the  velocity  away  from  the  leading  edge  gives 


A 


Uo 


a 


3  [Sa 

iV  27 


(2.6) 


The  invisid  how  around  the  parabola  is  conveniently  described  in  curvilinear 
coordinates.  If  the  transformation 


-  2r/)  , 

y  =  ^(h  +  1), 


(2.7) 


is  introduced,  then  the  surface  of  the  parabola  becomes  ?]  =  0  and  the  coordinate 
^  varies  from  — cx)  to  cx)  on  the  surface  of  the  parabola.  A  mapping  between  {x,y) 
and  (^,'/7)  is  shown  in  Figure  2.1.  The  surface  velocity  around  the  leading  edge  is 
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y 


Figure  2.1:  Parabolae  of  constant  ^  and  77,  with  the  airfoil  surface  defined  by  77  =  0 

=  (2,8) 

This  result  is  consistent  with  (2.4)  and  cte,  dehned  as  the  effective  angle  of  attack 
(see  comments  below  2. 21, and  also  2.41).  When  is  beyond  a  critical  value,  ap¬ 
proximately  1.16,  the  boundary  layer  separates  within  a  time  0(1),  scaled  using  the 
local  length  Vo  and  velocity  Uo-  The  streamwise  position  of  the  separation  point 
is  less  than  2  on  the  upper  surface  of  the  parabola.  This  is  relatively  close  to  the 
leading  edge  vertex  and  similar  results  are  expected  for  round-nosed  airfoils. 

The  vector  form  of  the  non-dimensional  Navier-Stokes  equations  for  incompress¬ 
ible  flow  is 


dv 

dt 


+  v  X  {V  X  v)  +  '^) 


-Vp  -  —V  X  (V  X  v). 
Re 


(2.9) 


V-v  =  0. 


(2.10) 


General  expressions  for  the  vector  operators  in  the  curvilinear  coordinates  (Ch^C) 
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are 


_  I  d  ^ 

V  =  T-^ei  +  —^62  +  7-^63, 

hi  Oq  h2  oTj  d(^ 

hiei  h2e2  he^ 

_  ^  -  I  <9  d  d 

V  X  V  = 


(9a;i  drj  d( 

hiU^  h2Urj  h^Ui^ 


(2.11) 


where  ei,  62  and  63  are  the  unit  vectors  in  each  of  the  coordinate  directions  and 
u^,  Uri  and  are  respective  velocity  components.  In  the  normal  direction,  suitable 
scaled  variables  are 

V  =  M^Re^/^  Y  =  T^Re^/^.  (2.12) 

For  the  two-dimensional  flow  held,  the  momentum  and  continuity  equations  can  be 
formulated  in  parabolic  coordinates  {^,'f])  as 


duf 


dt  +  i 


due  due 

— 


1  dp  ^  1  1  di^u^ 


dp  J  _|_  I  d^  Re  -I- 1  dp 

dp 


dp 


=  0 


d^  dp  + 1 

Let  be  written  as  u.  New  variables  are  introduced  through 


n  =  'Re^^‘^p\J -I-  1,  ^  ^  I  -I-  Idr  =  ^ 

Jo  2 


^\/?^n  +  sinh 
^nu 


3/2  ■ 


(2.13) 

(2.14) 

(2.15) 

(2.16) 
(2.17) 


^2  +  1 

V  2 

In  these  variables,  the  familiar  form  of  the  boundary  layer  equations  is  restored,  i.e. 


du  du  du 

dt  ds  ^  dn  ds  '  dri^  ’ 


dpoc  ^  d^u 


du  dv 
ds  dn 


Corresponding  boundary  conditions  are 

(■C  +  CXe 


U  ^  U cx^ 


+ 1)1/1 


as  n  ^  00 


u  =  V  =  0  at  n  =  0. 


(2.18) 

(2.19) 

(2.20) 
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(2.21) 


The  pressure  gradient  in  the  momentum  equation  (2.18)  is  given  by 
dpoo  _  dU^  Uoo  dUoc  _  {^  +  ae)il-aeO 

ds~dt  yen  ~  +  1)^/"  ’ 

The  external  velocity  appears  to  depend  only  on  the  angle  of  attack  «£.  This  angle, 
however,  is  not  merely  the  geometric  angle  of  attack,  but  an  effective  angle  of  attack 
including  effects  from  the  airfoil  geometry  and  unsteady  motion.  In  the  current 
approach,  when  the  leading  edge  is  approximated  by  a  parabola,  it  is  important  for 
the  inviscid  flow  calculation  to  evaluate  the  effective  angle  of  attack  ae,  see,  e.g. 
Section  2.3.2.  As  the  boundary  layer  nears  separation,  a  region  containing  the  zero- 
vorticity  line  thins  as  {tg  —  t)^^^  in  the  streamwise  direction  while  the  displacement 
thickness  grows  as  {tg  —  (Cowley  &  Van  Dommelen  1990);  tg  is  the  separation 

time.  Eulerian  calculations  encounter  great  difficulty  in  revolving  the  associated 
large  gradients  in  flow  variables.  Dynamically  rehned  grids  are  not  easy  to  implement 
since  there  is  no  clear  feature  in  the  Eulerian  held  that  guides  meshing.  A  Lagrangian 
approach  seems  to  be  the  optimal  adaptive  scheme,  and  the  variables  (n  and  v)  that 
become  singular  at  separation  do  not  appear  in  the  governing  equations.  For  various 
reasons,  Lagrangian  simulation  is  time  consuming.  It  is  convenient  to  switch  from 
an  Eulerian  simulation  to  a  Lagrangian  one  only  at  the  hnal  stage  of  the  simulation. 
In  the  Lagrangian  calculation,  huid  particles  are  labeled  by  new  coordinates 
At  the  initial  stage,  they  are  specihed  as  the  physical  positions  of  particles, 

sii:  V,  to)  =  i,  n{i,  ?7,  to)  =  ff  (2.22) 


Subsequently,  the  two  sets  of  coordinates  are  related  by  the  Jacobian  matrix. 


/(9s  (9s  \ 

di  dfj 
dn  dn 

\di  dp) 


(2.23) 


With  the  initial  conditions  (2.22),  the  continuity  equation  (2.19)  is  equivalent  to  the 
requirement  that  the  determinant  of  the  Jacobian  matrix  is  unity,  i.e. 


(9s  dn  (9s  dn 
df]  dfj 


(2.24) 


The  governing  equations  are  now 


du  (9  Poo  d‘^u  ds 

dt  ds  dn"^  ’  dt 

Derivatives  with  respect  to  n  can  be  represented  as 

(9  ds  d  ds  d 
dn  dfj  d^  dfj 


(2.25) 


(2.26) 
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which  is  then  substituted  into  the  momentum  equation  to  give 

du  (9 Poo  r  ds  d  ds  d  .2 

dt  (9s  ^  dfj  ^ 


(2.27) 


with 

|  =  (2.28) 

From  the  continuity  equation  (2.24),  the  normal  position  n  of  the  particles  is  found 
by  integration  employing  the  method  of  characteristics.  The  continuity  equation 
(2.24)  has  the  charicteristic  directions 


d.^  d?7  dn 

ds  ds  1  ‘ 

dfj  di 


(2.29) 


(.^,  rj)  are  interpreted  as  a  line  of  initial  locations  of  fluid  particles  arriving  at  constant 
s  at  a  hxed  time,  n  can  therefore  be  obtained  from 


'  wall 


s  -  +  s? 
s^  -r  s^ 


ds, 


(2.30) 


where  partial  derivatives  s^  and  Sf^  are  evaluated  at  each  time  step. 

A  singularity  occurs  when  a  stationary  point  appears  in  the  flow  held,  or  equiv¬ 
alently. 


ds  ^  ds 

=  0,  ,  -^  =  0. 


(2.31) 


df]  ‘  ‘  di 

As  there  are  no  singular  variables  in  (2.27),  Lagrangian  simulations  can  go  far  beyond 
the  separation  time  without  convergence  problems,  although  the  governing  equations 
are  no  longer  valid  because  of  the  initiation  of  strong  viscous/inviscid  interaction. 
Nevertheless,  Lagrangian  calculations  may  encounter  convergence  issues.  As  the 
particles  that  were  initially  close  to  each  other  migrate  apart,  their  physical  positions 
can  become  increasingly  distant  from  each  other,  and  more  iterations  are  needed  at 
successive  time  steps.  Because  of  this,  a  remesh  procedure  based  on  (2.24)  is  used 
when  necessary. The  singularity  that  the  Lagrangian  simulation  detects  may  be  a 
point  of  inhnite  displacement  when  the  boundary  layer  separates,  or  a  non-smooth 
behavior  on  a  streamline.  A  discussion  is  presented  in  Chapter  3. 


2.3  Boundary  Layer  Control 

2.3.1  External  Pressure  Field  Manipulation 

For  a  thin  airfoil  undertaking  an  unsteady  maneuver  (see  Figure  2.2),  it  is  con¬ 
venient  to  dehne  the  airfoil  surface  by  writing 

f  =  +  v{x,t)  +  g{x,t)  (2.32) 
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where  fi{x,  t)  is  the  thickness  or  ’’fairing”  distribution,  u{x,  t)  is  the  camber  function, 
and  g{x,t)  represents  the  unsteady  maneuver.  '+'  and  '  denote  upper  and  lower 
surfaces  respectively.  Temporal  changes  in  represent  a  flexible  surface.  Similarly, 
temporal  changes  in  u  can  be  used  to  describe  airfoils  with  drooped  leading  edges, 
and  even  flexible  surfaces.  The  time-dependent  form  g  can  allow  unsteady  maneu¬ 
vers.  In  what  follows,  Uoo  is  the  free  stream  speed,  c  is  the  chord  length,  e  is  the 
thickness/length  ratio,  a  =  0(e)  is  the  angle  of  attack,  and  'y{x,t)  is  the  vorticity 
strength  in  the  wake.  The  streamwise  velocity  for  unsteady  flow  around  the  airfoil 
can  be  expanded  as  (Zalutsky  2000) 

q=l  +  eqt  +  ...  (2.33) 


where  etif  is  the  hrst  order  perturbation,  given  by 


1 


d/U  d/U 

±__1  ll-x 

Try  ^-x  7r\ll  +  xjyi-^ 

-1  -1 


'1  +  e 


,dq  dq  dv  du. 

—  H - -  - ^ - 1 

^dt  dt  dC 


^-x 


i+t 


2tt\  1  -f  a; 


'e+l7(l  +  t-0 
e-1  e-a: 


(2.34) 


This  velocity  is  scaled  with  respect  to  Uoo  and  the  streamwise  coordinates  x  and  ^ 
with  respect  to  c.  Several  application  are  discussed  below. 
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2.3.2  Oscillating  Wings 


Oscillating  wings  are  often  used  to  study  dynamic  stall.  For  a  rigid  thin  airfoil 
rotating  around  its  center,  the  surface  function  is 

/  =  ±/i(a;)  +  u{x)  —  a{t)x  (2.35) 

The  unsteady  flow  velocity  on  the  surface  of  a  thin  airfoil  is  therefore, 

r  1 


g  =1  +  eq^  + 


1  +  6 


h'(0 


TT  i  ^+1 
-1 


± 


1 

±  — 


1  +  x  J  V  1-^ 
-1 


+  «(«)) 


?+i 


(2.36) 


■27rVl  +  a;y  y  I-  ^ 

-1 


e+1 


In  order  to  match  the  above  expression  with  the  leading  edge  velocity  (2.4),  a  new 
scaled  coordinate  and  scaled  time 

x={X  +  l)/ro,  T  =  — ,  with  ro  =  Ro(?,  (2.37) 

ro/Uoo 

are  introduced.  Here  Tq  is  the  nose  radius  and  Ro  =  64/27  for  the  Joukowski  airfoil 
(2.1).  The  dominant  behavior  of  the  surface  velocity  near  the  leading  edge  is. 


g  ~  1  ± 


RoX 


+  a{t)} 


di 


i+t 


t-i 


2 


+  -  /  7(l  +  t-0 


di 


Comparing  this  form  with  equation  (2.4)  shows  that, 
2 


V+^j’ 

(2.38) 


ctp  = 


TT 


1  l+t 

(  {d(«)«  +  o(i)}  +  l  f  1(1  + t-i) 


if, 


-1 


2  2 


(2.39) 


where  a^.  is  an  effective  angle  of  attack  and  a  is  a  geometric  angle  of  attack.  From 
the  conservation  of  global  circulation  (Zalutsky  2002) 


Fn  —  —2 


-1 


=  -2 


-1 


1±1, 

1-e 

^1+i 

i-e 


a{t)  +a{t)i)di 


i+t 


1 

i+t 


'ill 

e-i 


7(l  +  t-0  d^, 


(2.40) 


{a{t)  +  a{t)  d^  -  /  7(l  +  t-0 


d^ 


+  0(r, 
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where  Fq  is  the  global  circulation  around  the  wing  and  wake.  Combining  (2.40)  and 
(2.39)  leads  to 

«.(r)  =  {a  -  +0(4).  (2.41) 

For  the  situation  when  the  airfoil  oscillates  harmonically  around  a  constant  angle 
(Uo,  with 

a{T)  =  Rq  {tto  +  rotti  (1  -  coscoqT)}  ,  (2.42) 

the  effective  angle  of  attack  is  then 

ae{T)  =  tto  -  sincooT.  (2.43) 

Eulerian  and  Lagrangian  calculations  results  are  discussed  in  Chapter  3. 


2.3.3  Drooped  Leading  Edge 

A  drooped  leading  edge  is  rotated  about  the  quarter-chord  point  of  the  airfoil. 
The  geometric  angle  of  attack  is  a  and  the  drooping  angle  is  (5.  The  unsteady 
surface  function  is, 

+  iy{x) -a{t)x  + (3{t){x  + 1/2),  for  a;  < -1/2 
|^±/i(a;)  -I-  iy(x)  —  a{t)x,  for  x  >  —1/2 

The  horizontal  surface  velocity  component 


IS 


la.  la.  I  1 


-1 


± 


± 


1  +  X 


-1/2  , - 

y  Y  ^  ^  +  «(^)  +  /3^  +  Pii)) 


1  4-  a; 


-1/2 


(“^  +  “(c) 


± 


27r  V  1  +  a; 


'^+C(i+«-o 


-1 


(2.45) 


1-1"  ''1+1 

Following  the  procedure  sketched  in  (2.36)-(2.41),  the  effective  angle  of  attack 
be  obtained  as 


can 


hi  I  da  1  I 
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As  an  example,  let  a  be  fixed  at  ao,  and  let  the  drooped  nose  oscillate  around  a 
constant  angle  Pq.  The  oscillation  amplitude  is  (3i  and  frequency  is  Uq 

P{T)  =  rI  |/3o  +  roA  (1  -  coso^oT)}  .  (2.47) 

Thus 

,  ,  ,1  ^/3\  /I  S^/S/dlC^O  .  m  /rx 

ae{T)  =  cto  +  (x  +  75“)/^o - 3 - sinc^oT.  (2.48) 

O  Zn  OTT 

This  case  is  equivalent  to  the  oscillating  wing  result  (2.43),  but  with  a  higher  ao 
and  lower  di.  Similarly,  for  a  fixed  drooped  nose  P  =  Pq  with  the  airfoil  oscillating 
about  the  mid-chord  point,  from  (2.46) 

aPT)  =aQ  +  Po-  sincuoT.  (2.49) 

This  case  corresponds  to  a  higher  ao  than  (2.43).  Based  on  the  above  discussion, 
separation  on  a  maneuvering  airfoil  with  a  drooped  nose  can  be  studied  as  a  special 
case  of  an  oscillating  airfoil. 

2.3.4  Suet  ion/ Inject  ion 

Suppose  that  an  injection  slot  is  located  on  a  section  of  the  wall  from  C  to  D  and 
a  suction  slot  from  A  to  B.  This  situation  is  sketched  in  Figure  2.3.  The  velocity 


Figure  2.3:  Schematic  diagram  of  the  passive  device.  The  suction  slot  is  between  A 
and  B  while  the  injection  slot  is  between  C  and  D.  Both  the  injection  and  suction 
are  taken  here  to  be  normal  to  the  surface. 

distribution  in  the  suction  slot  is  taken  as 

Vs  =  sin’"  a)  ’ 
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while  the  velocity  distribution  in  the  injection  slot  is 

Vi  =  14i  sin”  .  (2.51) 

n  =  0  corresponds  to  uniform  suction  and  injection,  n  =  3,  referred  to  as  cubic 
suction  and  injection,  can  provide  smooth  influx  into  the  boundary  layer.  The  air 
flow  rate  in  the  suction  slot  should  equal  that  of  the  injection  slot.  Thus 

f  K;,sin”  j  f4,,sin”  (2.52) 

which  gives 

Vwi  =  Kos  (  £)  _  (^  )  ■  (2.53) 

Over  AB,  the  average  pressure  is  higher  than  over  CD;  flow  is  ejected  from  AB 
to  CD.  When  suction/injection  is  implemented  on  the  surface  of  the  parabola,  the 
governing  equations  (2.18)-(2.19)  do  not  change.  For  boundary  conditions  on  the 
surface,  the  normal  velocity  component  at  suction/injection  slots  is  given  by  (2.50) 
and  (2.51)  and  is  zero  elsewhere. 


16 


Chapter  3 

Two-dimensional  Numerical 
Results 


3.1  Oscillating  Wings 

Separation  on  oscillating  wings  is  simulated  at  various  frequencies,  amplitudes 
and  angles  of  attack.  Without  oscillation,  the  boundary  layer  separates  when  the 
scaled  angle  of  attack  is  beyond  a  critical  value  1.1556  (Werle  &  Davis  1972,  Ruban 
1981).  However,  the  effective  angle  of  attack  during  the  oscillation  involves  both 
steady  and  time  dependent  terms.  As  shown  in  (2.43) 

,rp\  dlC^O  . 

ae[T)  =  oo - —  smcuoT. 

Clearly  at  high  frequencies,  the  time  dependent  term  has  a  signihcant  effect.  From 
Tables  3. 1-3.4,  it  can  be  seen  that  the  separation  time  varies  considerably  for  os¬ 
cillating  wings.  Previously,  Zalutsky  and  Walker  (2003)  had  suggested  that  at 
cto  =  2.0,  di  =  1.0,  as  ojq  increases,  the  separation  time  first  fluctuates  around 
some  value  and  then  seems  to  blow  up  when  ujq  >  20.  Therefore  extensive  studies 
are  carried  out  for  various  angles  of  attack  with  oscillation  amplitudes  from  0.25 
to  4.0  and  frequencies  as  high  as  2000.  Very  small  steps  have  to  be  used  when  the 
frequency  becomes  high.  The  following  time  steps  are  used  in  the  simulation 


dt 


0.0005,  uo  <  50 

0.0002,  50  <  cuo  <  100 

0.00002,  100  <Uo<  1000 

0.00001,  1000  <  uo 


(3.1) 


At  large  frequencies,  it  may  take  several  million  steps  for  just  one  case  before  the 
simulation  fails  to  converge.  For  ckq  =  4.0,  the  proposed  case  studies  are  completed. 
The  results  are  plotted  in  Figure  3.1.  In  Figure  3.1(a),  for  di  =  1.5,  the  separation 
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time  T  increases  monotonically.  At  Uq  =  2000.0,  the  boundary  layer  separates 
at  T  =  8.46,  which  is  more  than  hve  times  the  non-oscillating  value.  In  Figures 
3.1(b)-(e),  the  separation  time  T  hrst  increases  with  ci^o  but  then  decays  drastically 
after  uq  exceeds  a  certain  value.  The  spike,  which  corresponds  to  the  maximum 
separation  time  T,  occurs  at  lower  frequencies  as  di  increases.  From  Table  3.1,  at 
tto  =  4.0,  di  =  2.0,  the  maximum  separation  time  T  =  20.68  when  ujq  =  700.  For 
Oil  =  1-5,  the  maximum  separation  time  T  =  90.09  with  ujq  =  90.  It  is  interesting  to 
note  that  T  is  very  sensitive  to  both  ujq  and  di.  After  T  reaches  a  maximum  value,  a 
small  increase  in  o^o  usually  causes  a  substantial  drop  in  T,  as  shown  in  Figure  3.1(e). 
The  benehcal  effect  of  oscillation  seems  to  occur  only  for  di  <  1.5.  Further  increases 
in  Oil  over  most  of  the  frequency  spectrum,  expediate  separation.  Spikes  can  also 
be  seen  in  Figures  3.1(f)-(h),for  di  >  1.5.  The  computational  results  indicate  that 
if  an  airfoil  oscillates  at  a  tuned  frequency  and  amplitude,  the  separation  time  can 
be  delayed  by  a  factor  of  more  than  hfty,  although  the  averaged  angle  of  attack  is 
high.  In  the  simulation,  ujq  is  scaled  with  respect  to  a  reference  value  1/to,  where 
(see  2.37) 

(3.2) 

and  its  typical  size  is  10“"^  for  a  rotating  blade.  In  order  for  the  oscillation  to  be 
effective,  the  actual  frequency  value  should  be  as  high  as  10^  ~  10®s“^,  which  is  too 
high  for  current  engineering  practice. 

The  streamline  patterns  near  the  leading  edge  at  ccq  =  4.0,  di  =  1.5,  ojq  =  2.0  and 
cco  =  4.0,  di  =  1.5,ci;o  =  50  are  shown  in  Figures  3.2  and  3.3.  By  comparing  these 
two  hgures,  there  are  clearly  two  different  mechanisms  that  lead  to  the  failure  of  the 
simulations.  In  Figure  3.2(a), the  low  frequency  case,  streamlines  bend  slightly  near 
the  leading  edge  at  time  t  =  1.  From  Figure  3.2(b),  it  can  be  seen  that  a  bubble, 
or  flow  recirculation,  forms  on  the  upper  surface  and  the  skewing  of  streamlines 
becomes  pronounced.  In  Figure  3.2(c),  the  bubble  is  stretched  along  a  line  across 
which  streamlines  become  discontineous.  The  boundary  layer  has  separated  and 
large  fluctuactions  in  the  streamline  patterns  are  observed.  For  this  low  frequency 
case,  separation  is  initiated  by  the  appearance  of  a  zero-vorticity  line  and  culminates 
in  the  vorticity  spike.  The  high  frequency  case  is  different.  In  Figure  3.3(a),  an 
elongated  bubble  can  be  seen  on  the  surface  at  t  =  1.  Due  to  the  motion  of  the 
oscillating  airfoil,  streamlines  are  no  longer  smooth,  but  small  amplitude  oscillations 
appear  on  streamlines  upstream  of  the  leading  edge.  The  bubble  disappears  in 
successive  steps,  as  seen  in  Figures  3.3  (b)-(d),  and  the  fluctuation  amplitudes  on 
the  streamline  patterns  decay.  At  t  =  10  in  Figure  3.3(e),  a  discontinuity  appears  on 
the  streamlines  and  its  position  should  correspond  to  the  maximum  and  minimun  of 
the  oscillation  angles.  In  Figure  3.3  (f),  when  t  =  12,  streamlines  are  substantially 
perturbed  by  the  airfoil  maneuvers.  In  Figure  3.3  (g),  a  large  oscillation  is  seen  in 
the  flow  patterns  and  the  boundary  layer  separates  at  t  =  13.97.  Notice  that  for  this 
case,  flow  reversal  is  no  longer  seen  after  t  =  1  and  the  simulation  diverges  because 
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the  flow  fleld  is  subject  to  very  strong  disturbances  over  a  long  time. 

Lagrangian  methods  were  used  for  both  low  frequency  and  high  frequency  cases. 
At  low  frequencies,  ojq  =  2.0,  both  Lagrangian  and  Eulerian  results  agree  and  give 
a  separation  time  around  T  =  2.6.  For  the  high  frequency  case  with  ujq  =  50.0, 
the  Lagrangian  calculations  initiated  at  different  starting  times  quickly  terminate. 
The  Eulerian  calculations  suggest  T  =  13.9.  It  is  conjectured  that  Lagrangian 
simulations  are  very  sensitive  to  high  frequency  irregularities  that  occur  on  the 
steamlines. 
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(a)  (0  (frequency) 


1  10  1  00  1000 


(b)  (0  (frequency) 

Figure  3.1:  Separation  times  for  an  oscillating  wing  at  various  frequencies  and  ampli¬ 
tudes.  Averaged  angle  of  attack  ao  =  4.0,  and  oscillation  amplitudes  (a)  di  =  0.25 
,  (b)  di  =  2.0 
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(c)  CO  (frequency) 


1  10  1  00  1000 


(d)  (0  (frequency) 

Figure  3.1:  Separation  times  for  an  oscillating  wing  at  various  frequencies  and  ampli¬ 
tudes.  Averaged  angle  of  attack  ao  =  4.0,  and  oscillation  amplitudes  (c)  di  =  0.75 
,  (d)  di  =  1.0 
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(e)  ®  (frequency) 
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Figure  3.1:  Separation  times  for  an  oscillating  wing  at  various  frequencies  and  am¬ 
plitudes.  Averaged  angle  of  attack  ao  =  4.0,  and  oscillation  amplitudes  (e)  di  =  1.5 
,  (f)  di  =  2.0 
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(g)  « (frequency) 


1  10  1  00  1000 


(h)  (0  (frequency) 

Figure  3.1:  Separation  times  for  an  oscillating  wing  at  various  frequencies  and  am¬ 
plitudes.  Averaged  angle  of  attack  ao  =  4.0,  and  oscillation  amplitudes  (g)  di  =  3.0 
,  (h)  di  =  4.0 


23 


Figure  3.2:  Temporal  development  of  the  instantaneous  streamlines  for  =  4.0, 
OL\  =  1.5,  and  uj  =  2.0  in  physical  coordinates  [x^y]  at  times  (a)  t=1.0,  (b)  t=2.0 
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Figure  3.2:  Temporal  development  of  the  instantaneous  streamlines  for  oq  =  4.0 
di  =  1.5,  and  oj  =  2.0  in  physical  coordinates  (x,  ?/)  at  the  time  (c)  t=2.739 


Figure  3.3:  Temporal  development  of  the  instantaneous  streamlines  for  ao  =  4.0 
di  =  1.5,  and  oj  =  50.0  in  physical  coordinates  {x,y)  at  the  time  (a)  t=1.0 
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Figure  3.3:  Temporal  development  of  the  instantaneous  streamlines  for  ao  =  4.0, 
di  =  1.5,  and  oj  =  50.0  in  physical  coordinates  {x,y)  at  the  time  (b)  t=2.0,  (c) 


t=3.0 
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Figure  3.3:  Temporal  development  of  the  instantaneous  streamlines  for  ao  =  4.0, 
di  =  1.5,  and  oj  =  50.0  in  physical  coordinates  {x,y)  at  the  time  (d)  t=6.0,  (e) 


t=10.0 
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Figure  3.3:  Temporal  development  of  the  instantaneous  streamlines  for  ao  =  4.0, 
di  =  1.5,  and  a;  =  50.0  in  physical  coordinates  {x,y)  at  the  time  (f)  t=12.0,  (g) 


t=13.97 
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Table  3.1:  Calculated  results  for  the  separation  time  T*  and  separation  location 
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Table  3.2:  Calculated  results  for  the  separation  time  T*  and  separation  location 
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Table  3.3:  Calculated  results  for  the  separation  time  T*  and  separation  location 
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0.5 


1 


2 


5 


10 


15 


20 


25 


30 


35 


40 


45 


50 


55 


60 


65 


70 


75 


80 


85 


90 


95 


100 


200 


300 


400 


500 


600 


700 


800 


900 


1000 


2000 


15.23 


15.47 


14.37 


15.32 


16.97 


16.44 


16.04 


16.40 


16.64 


16.21 


16.24 


16.28 


16.43 


16.55 


16.80 


16.92 


17.05 


17.15 


17.28 


17.43 


17.65 


17.98 


18.10 


18.15 


23.13 


34.39 


45.40 


16.71 


15.97 


14.08 


16.63 


20.30 


19.59 


19.83 


20.29 


21.21 


21.85 


22.71 


23.93 


25.38 


26.17 


28.27 


29.81 


31.55 


32.23 


35.42 


39.02 


40.53 


42.77 


47.34 


53.33 


18.07 


15.40 


14.23 


19.69 


22.73 


24.16 


26.17 


28.75 


33.04 


37.43 


43.03 


50.31 


73.94 


96.03 


120.3 


137.0 


148.2 


155.9 


160.2 


162.9 


164.3 


163.3 


161.6 


20.06 


14.94 


14.40 


22.97 


26.63 


31.65 


38.22 


49.20 


65.72 


96.83 


167.19 


212.96 


238.72 


247.52 


258.24 


254.9 


22.29 


13.52 


7.69 


4.19 


42.02 


72.54 


221.2 


446.5 


473.2 


483.4 


438.3 


404.7 


377.6 


360.9 


356.6 


0.805 


0.776 


0.620 


0.527 


0.433 


0.409 


0.154 


2.0 

3.0 

4.0 

21.97 

21.13 

20.48 

12.74 

11.68 

11.05 

6.914 

6.155 

5.742 

42.75 

3.212 

2.957 

74.95 

3.189 

1.240 

3.487 

0.975 

0.653 

1.944 

0.651 

0.430 

5.831 

0.375 

0.319 

4.912 

0.284 

0.257 

3.778 

0.235 

0.214 

1.091 

0.200 

0.183 

0.642 

0.175 

0.160 

1.546 

0.151 

0.143 

0.387 

0.138 

0.129 

0.239 

0.127 

0.115 

0.220 

0.117 

0.106 

0.295 

0.116 

0.098 

0.186 

0.097 

0.091 

0.173 

0.092 

0.085 

0.163 

0.087 

0.080 

0.154 

0.080 

0.075 

0.150 

0.077 

0.071 

0.208 

0.072 

0.067 

0.163 

0.068 

0.064 

0.066 

0.035 

0.031 

0.043 

0.023 

0.021 

0.032 

0.018 

0.016 

0.0260 

0.014 

0.013 

0.022 

0.012 

0.011 

0.019 

0.010 

0.009 

0.016 

0.009 

0.008 

0.014 

0.008 

0.007 

0.013 

0.007 

0.006 

0.006 

0.003 

0.003 

Table  3.4:  Calculated  results  for  the  separation  time  T*  and  separation  location 
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3.2  Suction/Injection 

Table  3.5  (a),(b)  show  calculated  separation  times  for  uniform  suction/injection. 
The  distances  of  the  slot  locations  from  the  leading  edge  vertex  are  measured  us¬ 
ing  (2.16).  From  the  tables,  it  can  be  seen  that  increasing  the  injection/suction 
strength  can  initially  delay  separation.  As  the  injection  strength  increases  further, 
the  separation  time  decreases  and  can  become  smaller  than  the  no-control  value  of 
5.8.  The  locations  of  the  suction/injection  slots  can  also  have  a  considerable  effect 
on  separation.  In  particular,  there  are  optimal  locations  for  leading  edge  suction 
slots  that  almost  double  the  no-control  value  of  the  separation  time  (see  Tables  3.5 
a,b). 


14 

0.1 

0.2 

0.5 

1.0 

CD 

AB 

(0.05,0.4) 

(0.5,1. 5) 

5.812 

5.867 

6.202 

6.643 

(0,0.3) 

(0.5,2) 

6.623 

7.071 

7.011 

7.871 

(0.1,0.4) 

(0.5,2) 

6.57 

7.069 

7.063 

7.875 

(0,0.4) 

(0.5,2) 

6.594 

7.091 

7.135 

7.803 

(-0.1,0.4) 

(0.5,2) 

6.627 

7.052 

7.009 

7.772 

(-0.2,0.4) 

(0.5,2) 

6.663 

7.032 

7.059 

7.967 

(0,0.4) 

(0.5,2.5) 

6.511 

7.672 

8.842 

9.213 

(-0.1,0.4) 

(0.5,2.5) 

6.567 

7.742 

8.744 

9.664 

(0.05,0.35) 

(0.5,2.5) 

6.512 

7.67 

8.847 

7.299 

(-0.2,0.4) 

(0.5,2.5) 

6.635 

7.821 

8.796 

9.549 

(0,0.3) 

(0.5,3) 

6.211 

6.705 

7.652 

5.213 

(0,0.4) 

(0.5,3) 

6.175 

6.708 

7.213 

5.39 

(-0.1,0.4) 

(0.5,3) 

6.218 

6.699 

7.797 

6.116 

(-0.25,0.4) 

(0.5,3) 

6.279 

6.827 

8.572 

6.89 

(0.05,0.35) 

(0.5,3.5) 

6.012 

6.122 

5.993 

4.383 

(0,0.4) 

(0.5,3.5) 

6.012 

6.121 

6.032 

4.711 

(-0.15,0.35) 

(0.5,3.5) 

6.159 

6.256 

6.317 

5.288 

(-0.3,0.4) 

(0.5,3.5) 

6.018 

6.321 

6.653 

5.826 

Table  3.5:  (a)  Calculated  separation  time  for  suction/injection  at  an  angle  of  attack 
tte  =  2.0 
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1.5 

2.0 

3.0 

5.0 

CD 

AB 

(0.05,0.4) 

(0.5, 1.5) 

7.229 

7.849 

8.741 

4.656 

(0,0.3) 

(0.5,2) 

9.155 

5.666 

4.282 

1.26 

(0.1,0.4) 

(0.5,2) 

5.428 

4.209 

3.713 

1.012 

(0,0.4) 

(0.5,2) 

8.938 

5.415 

4.035 

2.719 

(-0.1,0.4) 

(0.5,2) 

8.831 

9.755 

4.599 

3.959 

(-0. 2,0.4) 

(0.5,2) 

8.908 

9.888 

7.044 

4.411 

(0,0.4) 

(0.5,2.5) 

4.793 

4.144 

3.735 

1.022 

(-0.1,0.4) 

(0.5,2.5) 

6.159 

4.523 

4.0 

3.642 

(0.05,0.35) 

(0.5,2.5) 

4.536 

4.025 

2.526 

0.786 

(-0. 2,0.4) 

(0.5,2.5) 

11.148 

5.469 

3.868 

3.814 

(0,0.3) 

(0.5,3) 

4.162 

3.844 

1.326 

0.712 

(0,0.4) 

(0.5,3) 

4.126 

3.769 

3.18 

0.811 

(-0.1,0.4) 

(0.5,3) 

4.575 

4.055 

3.691 

1.26 

(-0.25,0.4) 

(0.5,3) 

5.137 

4.404 

3.868 

2.507 

(0.05,0.35) 

(0.5,3.5) 

3.786 

2.835 

0.81 

0.682 

(0,0.4) 

(0.5,3.5) 

3.869 

3.703 

1.224 

0.654 

(-0.15,0.35) 

(0.5,3.5) 

4.335 

3.942 

3.688 

1.134 

(-0. 3,0.4) 

(0.5,3.5) 

4.796 

4.178 

3.788 

3.603 

Table  3.5:  (b)  Calculated  separation  time  for  suction/injection  at  an  angle  of  attack 
tte  =  2.0 
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Chapter  4 


Three-Dimensional  Boundary 
Layer  Separation 

4.1  Introduction 

Two  different  helicopter  flight  conditions  are  considered  for  the  three-dimensional 
boundary  layer  separation  studies:  hovering  and  forward  flight.  During  a  hovering 
flight,  in  a  frame  hxed  with  the  rotating  blade,  the  oncoming  velocity 

Uo  =  Dr  (4.1) 

in  the  streamwise  direction  is  time  independent.  The  lifting  blade  at  incidence  can 
be  represented  by  a  suitable  surface  vorticity  distribution  or  bound  circulation.  Ac¬ 
cording  to  Kelvin’s  theorem,  variation  of  the  bound  circulation  leads  to  the  shedding 
of  an  inboard  vortex  sheet  from  the  trailing  edge.  At  the  blade  tip,  the  bound  cir¬ 
culation  rolls  up  into  a  tip  vortex  causing  considerable  downwash  in  this  region.  In 
industry  a  cut-off  parameter  is  often  used  near  the  blade  tip  when  calculating  the 
thrust  coefficient  (Leishman  2000)  because  of  the  low  effective  angle  of  attack  due  to 
the  influence  of  the  tip  vortex.  Unlike  trailing  wakes  from  a  fixed  wing,  which  move 
away  after  shedding,  trailing  wakes  (including  the  inboard  vortex  sheet  and  the  tip 
vortex)  for  a  hovering  rotorcraft  remain  beneath  the  blades  and  impose  a  much 
larger  effect.  The  downwash  induced  by  trailing  wakes  is  discussed  in  Appendix  A. 
Forward  flight  is  a  more  complex  situation.  As  the  helicopter  has  a  forward  speed, 
the  oncoming  velocity  with  respect  to  the  blade  is  different  at  various  azimuthal 
locations  (shown  in  Figure  4.1).  An  advancing  blade  with  a  90°  azimuthal  angle 
experiences  the  highest  oncoming  velocity  while  a  retreating  blade  with  a  270°  az¬ 
imuthal  angle  encounters  the  lowest  oncoming  velocity.  The  oncoming  velocity  can 
be  written  as 


Uo  =  +  /ism(4/o  -|-  fit),  (4.2) 

Wo  =  jucos(^o  +  fit)-  (4.3) 
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Advancing  side 
(t>  =  90 


Figure  4.1;  The  advancing  and  retreating  side  of  a  blade  during  the  forward  flight. 

Here  fi  is  the  advance  ratio  and  4/  =  4/o  +  flt  is  the  phase  angle.  In  order  to  generate 
a  nniform  lift  and  balance  the  moment  on  the  blade,  the  pitching  angle  of  the  blade 
has  to  be  varied  throngh 

6{r,  T)  =  6*0  +  6*icCos(4/)  +  6isSin{'^)  +  ...,  (4.4) 

which  contains  the  collective  pitch  6o  and  the  first  harmonics  of  the  Fourier  series.  6ic 
is  the  lateral  cyclic  pitch  and  is  the  longitndional  cyclic  pitch.  On  the  advancing 
side,  the  effective  angle  is  fairly  low.  On  the  retreating  side  the  effective  angle  of 
attack  is  mucher  higher  to  compensate  for  the  rednced  flow  speed  (Gorton  &  Hoad 
2002).  Unsteady  boundary  layer  separation  on  the  leading  edge,  initiated  by  the 
adverse  pressure  gradient  at  a  high  angle  of  attack,  leads  to  the  formation  of  a 
dynamic  stall  vortex.  When  the  stall  vortex  resides  on  the  blade,  a  signihcant  gain 
in  lift  is  observed.  The  stall  vortex  quickly  detaches  from  the  airfoil  surface  and  is 
transported  downsteam.  Conseqnently,  the  blade  experiences  a  drastic  lost  of  lift 
and  a  severe  increase  in  pitching  moment.  Althongh  the  flow  is  nnsteady  during 
forward  flight,  the  external  flow  near  the  leading  edge  is  qnasi-steady  on  the  local 
time  scale.  The  spanwise  velocity  has  a  stronger  inflnnce  on  the  development  of  the 
bonndary  layer  compared  with  hovering  flight,  bnt  the  dominant  factor  is  the  high 
effective  angle  of  attack  on  the  retreating  blade  causing  unsteady  separation.  Some 
detailed  calculations  are  planned. 
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Figure  4.2:  The  leading  edge  surface  of  a  rotating  blade,  s  is  along  the  tangential 
direction,  and  n  is  along  the  normal  direction. 


4.2  Eulerian  Formulation  for  the  Boundary  Layer 

In  the  laboratory  frame,  the  Navier-Stoke  equations  for  incompressible  flow  can 
be  written  as  (see  2.9  &  2.10) 

dv  1  1 

—  +  V  X  (V  X  v)  +  -V(v  ■  v)  =  -Vp  -  —V  X  (V  X  v),  (4.5) 

at  2  Ke 


V  ■  V  =  0.  (4.6) 

As  shown  in  Figure  4.2,  the  frame  is  attached  to  a  blade,  and  rotates  about  the  y 
axis  with  an  angular  speed  hi.  The  incident  free  stream  is  in  the  x-direction. 

The  boundary  layer  equations  around  the  rotating  blade  are 


du*  du*  du*  du* 

I  I  I 

dt*  dx*  dy*  dz* 


dp* 

+  2nw*  +  v 
ox* 


d‘^u* 

dy*2  ’ 


(4.7) 


dw*  ^dw*  ^:dw* 
dt*  dx*  dy* 


+  w* 


dw* 

dz* 


dp* 

-^-2nu*  +  u 

dz* 


d^w* 

dy*^' 


du*  dv*  dw* 
dx*  dy*  dz* 


(4.8) 

(4.9) 
(4.10) 


where  2Vtw*  and  2Vtu*  are  due  to  the  Coriolis  force.  For  the  boundary  layer,  the 
pressure  gradients  are  dehned  by 
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(4.11) 


dp* 

dx* 

dp* 

dz* 


dU^ 

dt 

dW, 

dt 


+  u: 


+  u: 


) 

.dW*^ 

O 


'  dz* 


2m¥* 


dW* 

+  w* - +  2Q[/*  . 

“  dz*  oo 


Since  the  flow  structure  near  the  leading  edge  is  required,  the  following  nondimen- 
sional  variables 


t* 

To/ria’ 


Re  = 


TnQa 


e  = 


(4,12) 


are  introduced.  Here  c  is  the  chord  length,  a  is  the  radius  of  the  blade,  is  the  nose 
radius  at  the  leading  edge,  and  Re  is  the  local  Reynolds  Number.  Consequently, 
suitable  dimensionless  velocity  components  and  cooordinates  are 


u 


u  = 


Vta 


V  = 


X 


X  = 


y  = 


y 


Via 


w 


w  = 


hla’ 


(4.13) 


2:  e 


= 


Similar  to  the  two-dimenisonal  boundary  layer  formulation,  curvilinear  coordinates 
{s,n,z)  are  introduced  as  shown  in  Figure  4.2.  {s,n,z)  are  connnected  with  the 
Cartesian  coordinates  {x,  y,  z)  through 


ds  =  n  =  +  iJie^/2^  (4.14) 

while 

X  =  ^  (f  -  -  2p)  ,y  =  (4.15) 

(4.15)  is  the  parabolic  transformation  used  in  the  previous  chapter.  The  dimension¬ 
less  form  of  the  governing  equations  can  therefore  be  represented  as 


du 


dw 


du  du  du  dUoo  ,  _  dU^  dU^o 

ds  dn  dz  dt  ds  dz 


dw  dw 


u- 


ds 


dn 


-  2e{W, 
dw  dWr 


w 


a/^2  1  dn^ 


3‘^u 

+  ^  +  0(5), 


ew 


dz 


dt 
f  2e{U, 


,  dW^  ^  dW^ 

+  Uoo^ —  +ehFc 


ds 


u] 


a/^2  1  dn^ 


dz 
d^w 

+  ^+0{S), 


du  dv 
ds  dn 


dw 

dz 


(4.16) 


(4.17) 


(4.18) 
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where  S  =  Max{e^ ,1/ Re).  As  in  the  two-dimensional  boundary  layer  formulation 
(2.18  &  2.19),  the  above  equations  are  established  in  a  coordinate  system  that  fits 
the  surface  of  the  leading  edge  with  a  parabolic  cylinder  as  shown  in  Figure  4.2.  In 
the  above  equations,  {u,v,w)  are  velocity  components  along  {s,n,z),  respectively. 

For  the  current  domain,  s  ranges  over  (— cxd,  cx))  and  n  ranges  over  (0,  cx))  It 
is  convenient  to  map  this  domain  to  a  finite  one  and  cluster  points  where  large 
gradients  in  the  flow  variables  occur.  The  following  transformation  is  used. 


1  fl 

s  =  hs{s)  =  kstan{7i{s  — n  =  hn{n)  =  kntan{7r—),  z  =  z.  (4.19) 


In  the  new  computational  domain,  the  governing  equations  are 


du 

^d^u  ^du  ^du  ^^du  „ 

(4.20) 

'm 

—  -Ro-2  ^  ^  Zi''  +-^1+^1, 

dn^  ds  dn  dz 

dw 

^d^w  ^dw  ^dw  ^^dw  „ 

(4.21) 

dt 

—  -Ro-2  ^  Zl''  +-^2  +  r2, 

dn^  ds  dn  dz 

1  du  1  dv  dw 

(4.22) 

h'g  ds  h'j^dn  ^  dz 

where  the  coefficients  in  equation  (4.20)-(4.22)  are 


R^(h'Jn))-\  S  = 


u 


K{sy 

T  =  —  v)hn{n),  H  =  —ew, 


Ki  =  +2ew 


iFo  =  —2eu 


A/eTT’ 


(4.23) 

(4.24) 

(4.25) 


^  dU^  ,  dU^  ,  dU^  ^ 

hi  —  — — h  UoQ—z - h  eWoQ—z - 2eIFoo 


dt 


ds 


dz 


ve+T’ 


(4.26) 


r2  —  — - 1"  Uoo — ;; - h  ehFoo^:; - h  2eUo, 


dt  ^  ds 
and  hyS),  h'^ifi)  are  dehned  as 

dhs{s) 


dz 


VFTT’ 


h'is)  = 


dhnifl)  knTl  2t  ^ 


,  = /Cs7rsec^(7r(.s - )),  h'(h)  =  — —— = - sec“(7r— l 

ds  ^  ^  2”'  ’  dn  2  ^2’ 


(4.27) 


For  an  impulsively  started  boundary  layer,  the  initial  velocity  corresponds  to  the 
Rayleigh  solution 

u  =  Uoceri{p),  w  =  IFooerf(p),  (4.28) 
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where  the  Rayleigh  variable  p 
tained  as 


Woo  can  be  obtained  by  substituting  (4.29)  into  (4.17)  to  give 


Woo  — 


(4.30) 


Generally,  the  contribution  from  the  trailing  wake  has  to  be  included.  As  Woo  is  not 
bounded  when  s  becomes  large,  the  parabolic  surface  used  in  the  simulation  does 
not  extend  inhnitely.  The  boundary  conditions  are 


(m,  v,w)=0  at  n  =  0, 

u^Uoo(s,z),  w^Woo(s,z)  asp^l.  (4.31) 


in  which  Uoo{,s,z)  and  Woo{s,z)  are  defined  by  the  external  inviscid  solution.  The 
Rayleigh  solution  is  also  implemented  at  upstream  and  downstream  inhnity 

u  17ooerf(p),  w  IRooerf(p),  as  s  — ±1.  (4.32) 

At  both  ends  in  the  spanwise  direction  (the  blade  root  and  tip)  it  is  assumed 
that 

d 

—  =  0,  at  i  =  0, 1  (4.33) 

oz 

which  reduces  the  boundary  layer  equations  (4.16)-(4.18)  to  a  quasi  two-dimensional 
form  at  the  root  and  tip.  A  similar  approach  can  be  found  in  (Atik  2002).  Reduced 
equations  at  the  spanwise  boundaries  are 


du  ^,d‘^u  ^,du  ^,du  ^  , 

a+  ~  ^  a-2  ^  a-  a'- 

ot  on^  os  on 

(4.34) 

dw  ^,d^w  ^,dw  ^,dw  ^  , 

(4.35) 

1  du  1  dv 

(4.36) 

h'g  ds  h'^  dn 

where  the  coefficients  in  equations  (4.34)-(4.35)  are 

0  7/ 

,  s' =  T'={h'^{n)  v)K{n^ 

),  (4.37) 

^  ,  dUoo  ^  dUoo  ^  ,  dWoo  ,  dWoo 

J-  1  —  o,  +  Ooo  r.  ■,  i  2  —  Oj.  +  o 

dt  ds  dt  ds 

(4.38) 
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At  the  early  stage  of  the  simulation,  it  is  convenient  to  introduce  the  Rayleigh 
variables 


n  Y 

^  2v^’  2Vi' 

(4.39) 

Also  in  the  computational  domain,  the  transformed  variable  p  is 
by 

dehned  implicitly 

p  =  hp{p)  =  /cptan(7r|). 

(4.40) 

As  t  becomes  0(1),  the  boundary  layer  solution  approaches  the  Blasius  limit 
and  the  effective  boundar  layer  shrinks  to  the  wall  in  Rayleigh  variables.  Thus  the 
dependent  variables  (s,  p,  z)  in  the  simulation  are  switched  back  to  conventional 
variables  {s,n,z)  at  a  time  t  =  td  ^  0(1).  The  coeffcients  kn  in  (4.19)  and  kp  in 

(4.39)  are  related  by 

kn  —  ‘^'s/^kp. 

(4.41) 

The  governing  equations  in  Rayleigh  variables  are 

du  ^du  ^du  „ 

(4.42) 

~  -2  ^  ft'-  ft ^ 

at  op^  os  op  oz 

dw  ^dw  ^dw  „ 

(4.43) 

ot  op^  os  op  oz 

and  all  coefficients  remain  the  same  except  T,  which  is  now 

r  =  (p  +  h'Y)  -  2^ftv)h,(p)nu). 

(4.44) 

Initial  and  boundary  conditions  are  implemented  in  (4.28)-(4.32). 


4.3  Lagrangian  Formulation  for  the  Boundary  Layer 

As  discussed  at  the  beginning  of  this  chapter,  the  Eulerian  method  is  incapable  of 
resolving  the  large  gradients  in  the  streamwise  direction  when  the  singularity  begins 
to  emerge.  Therefore  it  is  replaced  by  a  Lagrangian  method  at  an  appropriate  time, 
say  t  =  to,  when  the  boundary  flow  is  still  smooth  but  close  to  separation.  The 
solution  is  continued  numerically  until  separation  occurs. 

Lagrangian  coordinates  (C,f/,C)  are  dehned  as  the  fluid  particle  positions  at 
t  =  to-  The  physical  position  (s,  n,  z)  and  velocity  (u,  v,  w)  of  each  fluid  particle  are 
functions  of 

The  Lagrangian  and  Eulerian  variables  are  related  at  t  =  to  through; 

{s,n,z)  =  {u,w)  =  {uo{i,fi,C),Wo{i,fi,C)),  att  =  to 


(4.45) 

(4.46) 
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The  subsequent  evolution  of  the  velocity  components  can  be  determined  from  the 
Lagrangian  boundary  layer  equations 


du 

W 


dU^  ,  9f/oo  ,  dUo. 
■ - +  +  eWoo - 


dt 


ds 


2e(W^ 


w 


dz 

d^u 


dw 

W 


dW^  ,  dW^  ,  dW, 

- - +  Uoo^ —  +  ehhoo - 


dt 


ds 


+  2e(f/oo  —  u) 


+ 


dz 

d'^w 


a/^2  X  dn"^ 


The  Lagrangian  derivative  d/dt'  is 


d 


d  d  d  d 

—  +  u—  +v—  +  w—. 
dt  ds  dn  dz 


(4.47) 


(4.48) 


(4.49) 


For  convenience,  it  will  be  denoted  as  d/dt  in  the  Lagrangian  analysis.  Note  that 
d/dn  should  not  appear  in  the  governing  equations;  it  can  be  replaced  using  the 
Jacobian  transformation 


A 

dn 


d  d  d 

{x^Zf,  -  Xf^Z^)^  +  {x^z-^  -  x^z^)-^  +  {Xi^z^  -  X^Zf,)-^, 


(4.50) 


so  that  the  current  particle  position  (s,  z)  can  be  found  numerically  from 


ds 

m 


=  u, 


dz 

m 


=  ew 


(4.51) 


To  render  the  Lagrangian  computational  domain  hnite,  the  following  transfor¬ 
mation  is  introduced 


f  =  Catan(A/2),  f)  =  Cfetan(7rr)/2), 

s  =  Ca  tan(7rs/2),  n  =  Cf,tan(7rh/2). 

The  governing  equations  for  (p  =  u  or  w  are 


d(f) 

'm 


^  ^  p  ,  o  Al 

g^2  +  dti^  ^  de  didv  dtidC 


^  d‘^(b  dd)  dd)  dd> 

+  Qcr — ^ - h  Rf — ^  +  Rf}~^  +  ^<t>  “t"  ^01 

^^d^dC  dri  ^(9C 


(4.52) 

(4.53) 


(4.54) 
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and  equations  (4.51)  for  the  streamwise  and  spanwise  particle  positions  are 
as 


ds  2  2/  -  /  \ 

■  cos  {ns/2)  ■  u, 


=  ew. 


dt  nca  ’  dt 

The  functional  coefficients  in  equation  (4.54)  are  dehned  by 


p~~  = 

J-  rjf]  ^  ? 


Q..  =  2nAB,  =  2C1BC,  Q,.  =  2nAC, 


P|  =  Q 


Rff  =  ff 


R^  =  n 


+  B^  +  1  +  dl(P  -  H) 

df]  <9C  ' 

+  B^  +  ]+B{F-  H) 

dr]  (9C  ' 

A^-Z^PA^C^4\+CiF-H) 

dr]  dC  '  ^  ' 


=  +2ew 

=  —2eu 


v/FTT’ 


veTT’ 


)  =  u, 

=  w, 


r  _  ,  TT  o  u/ 

r«s  —  — ;r- — h  Uoo—F. - 1“  elToo^:; - 2efTo, 


dt 

dW^ 


ds 


dz 


,  dW^  ,  dW^  ,  ^ 

+  Uoo — - 1"  elToo — - h  2eU^ 


dt  '  ds  '  dz  ' 

Similarly,  the  coefficients  in  the  above  equations  are 


n  = 


2  cos(7r.^/2)  cos(7rr)/2)  1 
ncb  cos{ns/2)  j  ’ 


=  w, 


=  w. 


^  ds  dz  ds  dz  ^  ds  dz  ds  dz  ^  ds  dz 

d(df]  df]d('  didC  dCdi'  dr) 


ds  dz 

di  df]  ’ 


P  =  .tan(.5/2)(i^  +  P^+(7i), 

\  di  dr]  dcJ 
H  =  n  (Aid.\i{n^/2)  +  B  tan(7rr)/2)^  . 


recast 

(4.55) 

(4.56) 

(4.57) 

(4.58) 

(4.59) 

(4.60) 

(4.61) 

(4.62) 

(4.63) 

(4.64) 

(4.65) 

(4.66) 
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Initial  conditions  are  : 

s{i,V,C'to)  =  i;  =  fi-  z{i,fiX,to)  =  C  (4.67) 

u{L  V,  C,  to)  =  Uoit  fj,  C);  wX,  fj,  C,  to)  =  WoiL  V,  C),  (4-68) 

Boundary  conditions  at  downstream  and  upstream  infinity  are 

u  Uoo{s,  z)eii{p),  w  II4o(-s,  ;2)erf(p),  as  ±1  (4.69) 

Note  that  Uoo{x,  z)  and  Woo{x,  z)  have  to  be  updated  with  newly  computed  variables 
X  and  at  each  time  step.  Also 


M  =  0,  V  = 

0,  ,  tc  =  0 

at  r)  =  0, 

(4.70) 

u  =  U^{s,z), 

W  =  Woo{s,z) 

at  r)  =  1, 

(4.71) 

KLv,ct)  =  L 

z{i,v,Lt)  =  c 

at  r)  =  h  =  0, 

(4.72) 

s{i,v,Lt)  =  ±1; 

z{i,fi,Lt)  =  c 

as  i  ±1, 

(4.73) 

Infiow  and  outfiow  across  the  root  and  tip  boundaries  does  have  a  weak  effect 
on  the  interior  boundary  layer  separation.  However,  this  problem  is  expensive  to 
handle  in  the  Lagrangian  formulation.  For  convenience,  in  some  cases  the  spanwise 
velocity  is  assumed  to  be  zero  at  the  blade  root  and  tip.  Thus  at  the  root,  tc  =  0 
a.t  z  =  (  =  0,  while  at  the  tip,  w  =  0a,tz  =  (  =  l.  Using  this,  x  and  u  can  be 

d 

evaluated  at  planes  i  =  0, 1  with  the  assumption  —  =  0  (see  4.33).  Following  a 

dC 

procedure  similar  to  that  for  the  Eulerian  calculations,  the  governing  equations  at 
the  tip  and  root  are 


du 

'm 


d‘^u  d‘^u  d'^u 

p'^  +  g'4—  + 


„,du  ^,du  ^  , 

+  S'^+T'—  +  TX 

di  drj 


(4.74) 


and 


for  the  streamwise  particle 
by 


f)  c  2 

—  = - cos^(7rs/2)  ■  u,  (4-75) 

at  TlCa 

positions.  The  functional  coefficients  in  (4.74)  are  defined 


g' 


-2H' 


ds  ds 

di  dfj  ’ 


(4.76) 
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where 


S'  =n' 


-Q 


ds  1 

^  d^s 

dfj  ' 

yd^dfj 

ds  1 

^  d'^s 

to 

1 

^>ds 

/  d^s 

^  di 

\d^dfj 

f  d'^s 

dv 


ds 

dfj 


d  5 

7rtan(7r?7/2)—  )  , 


d  S 

7rtan(7r^/2)^ 

<9^ 


ds 


T.  ,  _  dUoo  ,  rr 

-1-7/,  -  ^  .  “r  ^  C30 


dt 


ds 


n'  = 


I  2  cos(7r.^/2)  cos(7r^/2)  1 
V  TTCfe  cos(7rs/2)  j 


Boundary  conditions  become 


(4.77) 


(4.78) 

(4.79) 

(4.80) 


M  =  0  atJ7  =  0,  u  =  Uoo  at  r)  =  1,  (4-81) 

u  =  Uocerf{h/{2y/t),  s  =  dzl,  (4.82) 

Ki,V,to)=i  (4.83) 

u{i,  V,  to)  =  uoii,  v);  x{i,  r),  to)  =  i-  (4.84) 


For  these  equations,  the  solution  procedure  can  be  found  in  Zalutsky  (2000)  &  Atik 

(2002) 

A  singularity  occurs  when  a  stationary  point  emerges  in  the  boundary  layer.  In 
another  word,  the  coefficients  of  partial  derivatives  in  (4.50)  have  to  vanish  at  the 
same  time: 


dsdz  ds  dz  dsdz  dsdz  ds  dz  ds  dz\  ^ 

d(dfi  dfjd^  didC  dCdi'  dt)  d^  df) ) 


(4.85) 


4.4  Numerical  Schemes 

In  the  Eulerian  formulation,  the  variables  u{s,  n,  z,  t)  and  w{s,  n,  z,  t)  are  solved 
using  a  time-marching  approach.  For  the  Crank-Nicolson  method,  each  term  is  aver¬ 
aged  over  the  current  and  previous  time  steps  so  that  second  order  accuracy  in  time 
is  achieved.  Central  differencing  is  used  in  discretizing  spatial  derivatives  of  second 
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order,  while  first  order  terms  are  handled  by  upwind-downwind  differencing.  In  this 
way,  the  coefficient  matrix  is  diagonally  dominant  and  the  spatial  resolution  is  also 
second  order.  Details  of  the  numerical  scheme  can  be  found  in  Peridier  &  Walker 
(1989).  Although  the  Lagrangian  formulation  has  different  variables: 

V)  C)  i))  'S(C)  V)  Cl  V:  C)  t)j  fhe  same  discretization  ideas  apply.  Matrix  iter¬ 

ations  implement  the  ADI  method,  which  sweeps  the  computational  domain  along 
grid  lines.  The  variables  on  the  current  sweeping  line  are  regarded  as  unkown  while 
neighbouring  points  are  assumed  to  be  known.  With  boundary  conditions  imposed 
on  both  ends  of  the  line,  unknown  variables  are  found  using  Thomas  algorithm. 
The  sweeping  direction  can  hrst  be  chosen  as  the  X  direction  so  that  every  grid  line 
along  the  X  direction  is  succesively  covered.  Then  the  direction  can  be  switched  to 
Y  or  Z.  As  updated  values  of  variables  can  be  used  directly  during  the  sweep,  the 
solution  usually  converges  within  2  or  3  iterations  at  each  time  step. 
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Chapter  5 


Three-dimensional  Numerical 
Results 

5.1  Eulerian  Calculations 

5.1.1  External  Velocity 

External  velocity  components  in  the  computational  domain  are  plotted  in  Figure 
5.1.  For  an  effective  angle  of  attack,  a  linear  twist 

ae  =  4-3z  (5.1) 

is  assumed.  Downwash  from  the  trailing  wake  is  not  included.  The  ratio  e  =  ^  is 
chosen  as  0.001  and  the  chord  thickness/ length  ratio  is  assumed  to  be  0.1. 

In  Figure  5.1(a),  at  a  hxed  section  along  the  span,  Erst  increases  in  the  pos¬ 
itive  streamwise  direction  (from  the  lower  surface  to  the  upper  surface).  Beyond  a 
certain  maximum  value,  Coo  then  decreases  due  to  the  adverse  pressure  gradient.  In 
Figure  5.1(b),  Woo  has  a  very  slow  variation  over  most  of  the  domain  in  the  stream- 
wise  direction.  The  rapid  variation  near  s  =  0  or  1  is  due  to  the  transformation 
(4.36)  used.  The  actual  gradient  of  Woo  in  the  streamwise  direction  near  s  =  0  or  1 
is  small. 

5.1.2  Separation  Time 

Boundary  layer  separation  on  the  leading  edge  of  a  rotating  blade  is  weakly 
three-dimensional  due  to  the  large  spatial  disparity  between  the  streamwise  and 
spanwise  directions.  In  the  streamwise  direction,  the  velocity  component  has  a 
variation  0(1)  within  a  leading  edge  nose  radius.  Conversely,  the  spanwise  velocity 
component  has  a  small  magnitude  and  a  small  variation  over  the  length  of  the  blade. 
In  the  streamwise  momentum  equation,  to  leading  order,  the  contribution  from 
the  spanwise  direction  is  insignificant.  Strong  gradients  in  streamwise  and  normal 
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Figure  5.1:  (a)  The  streamwise  velocity  in  the  outer  flow. 


Figure  5.1:  (b)  The  spanwise  velocity  in  the  outer  flow. 
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velocity  components  are  developed  in  a  way  similar  to  two-dimensional  separation. 
It  is  expected  that  for  certain  simple  conditions,  such  as  hovering  flight,  separation 
location  and  time  can  be  predicted  from  the  two-dimensional  separation  results.  The 
relation  between  separation  time  and  the  effective  angle  of  attack  is  shown  in  Figure 
5.2.  Dimensionless  variables  used  in  the  two-dimensional  boundary  layer  studies  are 
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Figure  5.2:  The  relation  between  separation  time  and  the  effective  angle  of  attack 
for  two-dimensional  leading  edge  separation. 


2.5  3  3.5 

the  effective  angle  of  attack 


The  scaled  oncoming  velocity  is  always  one.  For  a  rotating  blade,  the  upstream 
flow  speed  varies  linearly  along  the  span.  As  a  preliminary  three-dimensional  case, 
suppose  the  effective  angle  of  attack  varies  along  the  span  according  to  a  linear  twist 
(5.1),  where  downwash  from  the  trailing  wake  is  not  included.  In  hovering  flight, 
the  upstream  flow  speed  is  linear 


Doo  =  Dr.  (5.3) 

At  each  cross-section,  a  local  separation  time,  scaled  with  the  local  upstream  velocity 
in  (5.3),  can  be  calculated  with  respect  to  the  effective  angle  of  attack.  This  trial 
case  is  displayed  in  Figure  5.3.  The  minimum  separation  time  predicted  based  on 
the  two-dimensional  separation  results  is  around  7.0  in  the  interval  G  [0.4,  0.5] 
along  the  span.  This  three-dimensional  separation  occurs  at  T  =  6.69. 
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Z  spanwise  coordinate 


Figure  5.3:  The  relation  between  separation  time  and  the  effective  angle  of  attack 
for  two-dimensional  leading  edge  separation. 

5.1.3  Displacement  Thickness 

Three  dimensional  Eulerian  calculations  have  been  carried  out  for  the  boundary 
layer  structure  valid  near  the  leading  edge  of  a  blade.  The  evolution  of  displacement 
thickness  for  the  effective  angle  of  attack  (5.1)  is  displayed  in  Figures  5.4(a)-(d).  In 
this  case,  the  Eulerian  simulation  fails  to  converge  at  t  =  6.679.  As  the  dimen¬ 
sionless  time  t  (scaled  with  respect  to  the  nose  radius  and  the  blade  tip  speed) 
approaches  this  value,  spiky  behavior  is  seen,  and  indicates  a  local  singularity.  The 
corresponding  spatial  location  is  taken  to  dehne,  approximately,  the  leading  edge 
separation  location.  Changing  the  computational  time  step  from  0.005  to  0.001  does 
not  influence  the  onset  and  location  of  the  singular  eruption. 

5.1.4  Zero  Vorticity  Iso-surfaces 

Three-dimensional  separation  theory  predicts  that  both  vorticity  components 
dvj  dw 

00 z  ~  —7:—  and  oo^  ~  t:—  vanish  when  separation  occurs.  For  a  rotating  blade  with 
oy  ay 

a  high  aspect  ratio,  however,  the  magnitude  and  direction  of  the  spanwise  velocity 
component  only  have  a  very  weak  influence  on  the  separation  location  and  time, 
but  they  do  play  a  role  in  dehning  the  zero  iso-surface  of  the  streamwise  vorticity 
component.  If  the  spanwise  velocity  flow  is  from  the  outboard  to  inboard  (or  cUa,  <  0 
initially),  then  both  zero  iso-surfaces  of  oox  and  ooz  appear  and  intersect  when  the 
boundary  layer  nears  separation,  see  Figure  5.5(a).  Otherwise,  if  the  spanwise 
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Figure  5.4:  Displacement  thickness  at  t=6.1,  6.3,  6.4,  6.5,  in  (a),(b),(c),  (d)  respec¬ 
tively. 
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Figure  5.5:  (a)  Zero  vorticity  iso-surfaces,  the  red  one  represents  and  the  yellow 
one  represents 

flow  is  from  the  inboard  to  outboard,  only  the  zero  iso-surface  of  ojz  is  found,  and 
the  intersection  between  =  0  and  0;^,  =  0  could  not  be  determined,  see  Figure 
5.5(b).  The  requirement  =  0  and  cUa,  =  0  for  a  rigorous  separation  criteria  comes 
from  the  idea  that  the  displacement  thickness  tends  to  inhnity,  or  in  the  continuity 
equation  y  ^  cx)  at  the  singular  point.  For  a  rotating  blade  with  a  large  aspect 
ratio  the  spanwise  velocity  has  a  small  contribution  to  the  continuity  equation.  An 
accurate  determination  of  zero  vorticity  surface  for  is  difficult.  Nevertheless,  the 
displacement  thickness  can  still  become  very  large  at  certain  points  and  leads  to 
failure  of  the  simulation. 


5.2  Separation  Control 

Uniform  suction  is  implemented  on  the  leading  edge  along  the  whole  blade  length. 
From  the  preliminary  results,  suction  can  considerably  delay  boundary  layer  sepa¬ 
ration  on  a  rotating  blade.  In  Table  5.1,  suction  starts  from  the  leading  edge  vertex 
and  extends  over  a  part  of  the  upper  surface.  E.g.,  0.002  in  Table  5.1  means  0.2 
percent  of  chord  is  subject  to  suction. 

In  general,  for  three-dimensional  flow,  passive  suction/injection  devices  can  still 
be  used,  but  injection  should  occur  only  where  the  flow  will  be  weakly  perturbed. 
This  gives  more  flexibility  to  the  control  procedure. 
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Figure  5.5:  (b)  Zero  vorticity  iso-surface,  the  yellow  one  represents  ujz- 


0.1 

0.2 

0.5 

1.0 

2.0 

3.0 

5.0 

10.0 

slots 

0.002 

6.745 

6.74 

6.83 

6.99 

7.295 

7.505 

7.875 

8.15 

0.01 

7 

6.7 

6.755 

7.735 

9.105 

10.175 

11.6 

8.41 

0.023 

7.145 

7.74 

9.31 

11.065 

16.115 

19.455 

22.51 

8 

0.040 

7.14 

7.76 

10.635 

17.925 

29.245 

37.86 

28.72 

8.085 

0.06 

7.14 

7.76 

10.68 

35.15 

>  60 

>  60 

28.2 

8.6 

0.088 

7.14 

7.76 

10.675 

57.725 

>  60 

>  60 

27.975 

9.615 

Table  5.1:  Calculated  separation  time  Tg  with  suction  I4,  applied  on  leading  edge 
of  a  rotating  blade 
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Appendix  A 


Lifting-line  Theory  for  the  Downwash  of  a 
Multi-bladed  Rotor  with  Twist  in  Hover 


Following  the  approach  used  by  Li  (2000),  the  downwash  distribution  along  hov¬ 
ering  blades  can  be  obtained  from  lifting-line  theory.  The  geometric  angle  of  attack 
a  =  a{r),  where  r  is  the  distance  along  the  blade;  r  =  a  is  the  blade  radius.  For  a 
rectangular  blade  the  bound  circulation  for  a  rectangular  blade  is  given  by 


F  =  7rc(i}ar  —  Vz), 


(A.l) 


where  c  is  the  constant  chord  length.  In  (A.l),  Vz  is  the  downwash  induced  by  the 
trailing  vortices.  The  helical  vortex  shed  between  a  section  (r',  r'  +  dr')  from  a  blade 
is  assumed  to  be  a  semi-inhnite  cylinder,  and  the  strength  of  the  vortex  cylinder  is 
obtained  by  averaging  the  helical  vortex  over  the  pitch  between  neighbouring  spirals. 
A  semi-inhnite  cylinder  of  radius  r',  composed  of  vortex  rings  of  constant  circulation 
per  unit  length  70,  induces  a  vertical  velocity 

{^,  for  r  <  r' 

for  r  =  r'  ,  (A. 2) 

0,  for  r  >  r' 


For  a  multi-bladed  rotor  there  are  n  blades.  The  strength  of  the  vortex  cylinder 
is  the  sum  of  the  individual  contribution  from  each  blade.  Notice  in  (A. 2),  a  vor¬ 
tex  cylinder  does  not  induce  velocity  at  outside  points.  Thus  contributions  to  the 
downwash  at  radius  r  in  the  rotor  plane  are  associated  only  with  cylinders  of  larger 
radius.  The  net  downwash  is 

Vz{r)  =  \  lei{r')dr'  +  (A.3) 

where  the  subscripts  i  and  t  represent  the  inboard  vortex  sheet  and  the  tip  vortex 
respectively.  In  one  rotor  revolution,  the  vortex  at  radius  r  is  swept  downward 
through  a  distance  2ttVz/^.  As  F  is  the  bound  circulation  along  the  blade,  the 

strength  of  the  shed  circulation  per  unit  length  along  the  rotor  is  — — .  Therefore, 

dr 

for  a  vortex  cylinder  with  a  radius  r,  the  azimuthal  component  for  n  blades  is 


le{r) 


nVt  dr 
2tiv z{r)  dr 


nVt 

2'KVz{r) 


da 

a  +  r- - 

dr 


1  dvz 
n  dr 


(A.4) 
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The  vertical  displacement  in  one  revolution  of  the  tip  vortex  can  deduced  from  (A. 2) 
as 


Az 


7^^ 

4  n  ’ 


(A.5) 


where  'jet  is  the  azimuthal  component  of  vorticity  at  the  tip.  Also  'yet  satishes 


,  ,  nT  2nTtQ 

Mr)  =  —  =  - 

A;^  Ti'yet 


which  gives 


Mr)  = 


nTt^l 
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Substituting  (A. 4)  and  (A. 7)  into  (A. 3)  leads  to 


(A.6) 


(A.7) 


nVPc 


vAr]  = 


,da  1  dv:, 


a 


dr'  n  dr 


dr' 

Vz{r') 


+ 


nTthl 

27r 


(A.8) 


Let  tto  be  the  reference  angle  of  attack,  assumed  to  be  of  the  same  order  as  the 
thickness/length  ratio  of  the  chord,  see  e.g.  (2.6).  After  introducing  the  following 
dimensionless  variables 


=  g  a  =  ii,  w(E)  =  M 

a  ao  izatto 


A  = 


aoa 


r  = 


nc  Qac 

the  dimensionless  downwash  can  be  written  as 


(A.9) 


At  the  tip. 


^  ^,d(a  dW(R') 

cy  lb 


dR  dR 


W{R)  V  2A 


dR  /  1 

^  +  (A-10) 


fT(l)  = 


+  8Ad(l)  -  1 


4A 


Differentiating  both  sides  of  (A.  10)  gives 


T^wdWiR)  ^da 


The  solution  of  the  above  equation  is 


W{R)  = 


VI  +  8AF  -  1 


4A 


where 


F  = 


da 

/ 

a  +  R 

Ir 

dR' 

dR'  +  «(!). 


(A.ll) 

(A.12) 

(A.13) 

(A.14) 
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If  the  rotor  is  untwisted,  or  a  =  /3o,  where  (3q  is  a  constant,  then 


W{R) 


8\I3qR  —  1 
4A 


If  the  rotor  has  an  ideal  twist,  or  a  =  /3o/-R,  then 


(A,15) 


W{R) 


\/\  -\-  8A/3o  —  1 

4A 


(A.16) 


which  is  consistent  with  the  discussion  in  Leishman  (2000).  If  the  rotor  has  a  linear 
twist,  oi  a  =  Pq  —  PiR,  then 


W{R)  = 

For  the  effective  angle  of  attack 


v/l  +  8A(/3oi?-/5ii?2)-i 


4A 


(Up  =  a  — 


VLr 


or 


die  =  ttO  /3o  -  A-R 


^l  +  ^\{(3oR-  (3iR^)-l 


A\R 


(A.17) 


(A.18) 


(A.19) 


in  the  linear  case.  The  derivative  of  cte  along  the  spanwise  direction  is 


^  ^  ^  v/l  +  8A(/3oi?-/3ii?2)-l  _  /do  -  2(3iR  ^ 

dR  4Ai?2  R^l  +  S\{j3oR- l3iR^)  J 

(A.20) 

At  the  blade  root  R  =  0,  ag  vanishes  and  its  derivative  along  the  spanwise  direction 
is  hnite. 
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